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Abstract. We analyze the Lattice Boltzmann method for the simulation of fluctuating 
hydrodynamics by Adhikari et ah [Europhys. Lett. 71, 473 (2005)] and find that it 
shows excellent agreement with theory even for small wavelengths as long as a sta- 
tionary system is considered. This is in contrast to other finite difference and older 
lattice Boltzmann implementations that show convergence only in the limit of large 
wavelengths. In particular cross correlators vanish to less than 0.5%. For larger mean 
velocities, however, Galilean invariance violations manifest themselves through errors 
of a magnitude similar to those of the earlier implementations. 



1 Introduction 

Fluctuations are important for many hydrodynamic phenomena, from colloid diffusion 
to phase-separation close to the critical point. Particle based methods such as Stochastic 
Rotation Dynamics |2j, Lattice Gas [3 J or Molecular Dynamics simulations |4] naturally 
give rise to stochastic noise. In contrast the lattice Boltzmann (LB), or finite difference dis- 
cretization of the Navier Stokes equations require fluctuations that have to be included 
manually. The guiding principle for doing this is the theory of the fluctuating Navier 
Stokes equations [5J. Despite the success of applying the Navier Stokes equations to 
very small-scale flows formally the hydrodynamic limit requires large wavelengths. For 
fluctuating hydrodynamics the constraint of large wavelengths becomes important and 
standard discretization will give results that are not in agreement with statistical physics 
for shorter wavelengths. For a detailed analysis of simulating fluctuating hydrodynam- 
ics using finite difference methods and some remedies to improve this situation see the 
recent manuscript of A. Donev 0. Similar deficiencies are found for implementations of 
fluctuating Navier Stokes equations using the Lattice Boltzmann approach introduced by 
Ladd 1 7]. It is, however possible to use a more fundamental approach to include fluctua- 
tions in the LB method. Adhikari et al. [TJ introduced noise on all nonconserved modes, 
not only the hydrodynamic ones, leading to a scheme which shows good agreement with 
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theory even for large wavelengths. Duenweg et al. rederived this noise implementation 
from detailed balance considerations of lattice gases |8|. Both approaches are numeri- 
cally identical. In this paper we study the degree of improvement achieved and show 
that many of the deficiencies that plague finite difference discretizations of fluctuating 
Navier Stokes equations are absent in this Lattice Boltzmann implementation as long as 
we consider a system with vanishing mean velocity. For large mean velocities Galilean 
invariance is violated and errors of a similar magnitude to the earlier implementations 
are observed. 



2 Fluctuating Lattice Boltzmann with Ghost Noise 



Following the derivation of Adhikari et al. [H we start with the Lattice Boltzmann equa- 
tion (LBE) 

f l (x+v l/ t+l)=f i (x / t)+J2A ij \f j (x / t)-ff(x,t)}+^(x,t). (2.1) 



Here the the /, are the particle densities at position x, time t associated with with veloc- 
ity Vj. A;/ is the collision matrix and are the noise terms. We use the standard local 
equilibrium distribution given by 



f?=p Wi 



1 + ? U - Vi + 2? (U - Vi)2 -2? U - U 



(2.2) 



which is the discretized version of a Maxwell distribution f^OH- In equilibrium the /, 
will fluctuate around this distribution. The noise terms must be chosen such that, 
in the case of isothermal Lattice Boltzmann (LB), the density p = YLifi an d momentum 
pu = YLifi^i are conserved, i.e. Yd^i = and Y^i^i v i = 0- Furthermore a proper fluctuation 
dissipation theorem (FDT) corresponding to the collision operator A« is obeyed. This 
implies that the are correlated. We can find a representation in which the noise terms 
are uncorrelated by transforming the LBE into moment space. The moments are given by 



M fl (x,t) = E<f;M 



(2.3) 



So far this is a standard Multi-Relaxation-Time (MRT) representation |ITT1 - [T3| . The back 
transform is given by fi(x,t)=Y^ a n"M a (x,t). However, in order to construct a proper FDT 
these transforms cannot be orthogonal as in other MRT methods I11IH2J , so here we have 
rf-^m a -. Instead the transforms are chosen such that 



'^Wim a i m h i = ymftij = 5' 



ab 



(2.4) 



with n" = Wjm" while maintaining a diagonal moment space representation of the col- 
lision operator A,y = — X^X^w^^wj'. Now the moment transformation matrices are 
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orthogonal with respect to the Hermite norm. Such transforms with weighted norms 
were proposed before HT414T6H in different contexts. The necessity of the Hermite norm 
is briefly outlined after Eq. (|2.7D below and allows for a convenient definition of the mo- 
ment space noise terms as independent random variables. We can now rewrite the 
collision term of the Lattice Boltzmann equation in terms of the moments M" as 

/,(x+v i; f+l) = |M fl ( X/ - 1 [M a (x,f) - M a '°(x^)] +?} • (2-5) 

Adhikari et al. (H then obtain the FDT by performing a Fourier transform of the fluctua- 
tions from the mean of the moments SM"=M a — (M fl ) . They then use the fc-independence 
of these for an ideal gas to obtain 

<£T ) = ^Xl^ 1 (SM a SM c ) . (2.6) 

One particular result of the derivation is that the moment fluctuations £ fl decouple be- 
cause 

(sM*6M b ) = EL m H( s f' s fj) 

i i 
' ;' 

= Y^ m "i m} iP w i 

i 

= p5 ab . (2.7) 

Here we used {5fi5fj) = fiSu with Sfj = fi—fi where /; is the spatially uniform global 
equilibrium distribution function fT7\ . Adhikari also assumed that u <C 1 so that /, = pzo,. 
This allows us to use the orthogonality relation of Eq. (12.41) in the last step of the calcula- 
tion above. For a different transformation we would obtain non-diagonal elements in the 
fluctuation matrix which will then require correlated noise terms which are more cum- 
bersome to implement. For practical applications it is important to note that the u <C 1 
condition for the noise introduces a non-Galilean invariant contribution. We comment on 
this in our validation section. Inserting Eq. (|2.7[) into Eq. (|2.6|) leads to a noise expression 
of 

? = ± i yJp(2T«-l)N, (2.8) 

where N is a random variable with zero mean and a variance of one. 

Note that the moments M a are chosen to include the hydrodynamic moments. In the 
isothermal case discussed here they are comprised of the conserved quantities p and ], 
and the stress modes II. The remaining degrees of freedom are often called ghost modes 
as they do not appear in the isothermal Navier Stokes equations. However, the key result 
of the Adhikari et al. paper [l) was that they need to be taken into account when including 
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noise. Thus we add noise on all non conserved quantities, i.e. stress and ghost modes, in 
Eq. OJ| according to Eq. Oil . 

In practice we implement this algorithm by calculating the moments by means of 
Eq. A2.3D , performing the collision on the moments, adding the noise term and then trans- 
forming back into /-space as indicated in Eq. (|2.5|) . The streaming step is then done in 
/-space. This algorithm is almost as efficient as the standard LB implementation. The 
additional computational cost for calculating the ghost modes and the random numbers 
results in a computational overhead of less than 20%. 



3 Correlators in a D2Q9 Implementation 



To evaluate this method we present results for the D2Q9 (two dimensions, 9 base ve- 
locity vectors) LB model. The results are similar for other models, in particular we also 
tested D1Q3 and D3Q15. As D2Q9 base velocity set we use {v t } = {(0,0), (1,0), (0,1), 
(-1,0),(0,-1), (1,1), (-L1), (-1,-1),(1,-1)} and the {«;,•} = {4/9, 1/9, 1/9, 1/9, 1/36, 
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Figure 2: S^(u x ) averaged over 2xl0 7 iterations in a T a = l for all a, V = 20 2 , fluctuating D2Q9 ideal gas 
without and with active ghost noise. Note that different scales are used to visualize the slight deviations seen 
in the ghost noise case. 



1 / 36, 1 / 36, 1/36} . The matrix elements m\ in transform (|2.3|) are then given by 
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(3-1) 



The corresponding elements of the back transform are defined by the requirement 
ft" = mfwi. The zeroth moment then is the density p, the first and second are (up to a 
factor) the components of the momentum, the third and fourth the components of the 
shear stress and the fifth resembles the bulk stress [18] . The remaining three moments 



6 



are the ghost modes. Thus the equilibrium moments M a ' = Yd m 1fi are 

M°'° = p 

M 1 ' = V3pu x 

M 2 ' = V3pUy 

M 3 - = \p{u\-u^ (3.2) 

M 4 ' = 3pw jMy 

M 5 ' = fp(u*+«J) 

M 6,0 = M 7 / = M 8,0 = _ 

We present here results for ^-independence of the moment fluctuations predicted by 
Eq. (12.711 . In particular we consider the normalized static structure factor 

S k (M a ) = N a {SM a (k)SM a (-k)) (3.3) 

where 5M a (k) =X^ x ^M fl (x)e~' kx is the discrete spatial Fourier transform of 5M a and £ x is 
understood to be the summation over all discrete lattice sites. The normalization constant 
N a such that Sk(M fl ) = 1 is equivalent to p. I. e. for the density N p = and velocity 

components N Ux = p V \ bT where k^T = | for the isothermal D2Q9 model employed. A 
value of 1 throughout fc-space for the structure factor of any of the moments given in 
Eq. I|3.2|) thus indicates agreement with Eq. I|2.7|) . The volume V is just the number of 
lattice points V = £ x 1 and the division by it is just a normalization artifact of the Fourier 
transform. 

According to the argument put forth in [ 1 ] we expect the mean square fluctuations 
of all moments M a to be unity throughout fc-space. For the density p this is confirmed 
to three orders of magnitude in Fig.[ljb) for S k (p) and in Fig.[2pb) for Sk(w. T )- We find 
similar agreement for all nine moments of the D2Q9 model. For comparison we set the 
noise on the non-hydrodynamic modes (M 6 ,M 7 ,M 8 ) to zero, recovering the original Ladd 
method 17] and, as seen in Figures[T{a) and|2{a), the lack of noise on the ghost terms leads 
to drastic deficiencies for large k x ,k y values. Note that there are no deficiencies in Fig. [T] 
for k x = and k y = 0. The reason is that the projection of the D2Q9 model onto one 
coordinate axis yields a D1Q3 model. The isothermal ideal gas D1Q3 model, however, 
only has one stress mode and no ghost modes and thus there is no difference between 
the Ladd and Adhikari implementations in these projections. This is again observed in 
Fig.[2£a) where Sk(«x) exhibits white noise along the k x axis even in the absence of ghost 
noise. Motivated by private communication with A. Donev who is developing a general 
finite volume scheme to solve the fluctuating Navier Stokes Equations |6J based on a 
third order Runge-Kutta integrator we also measured the cross correlator 

R*{u x ,u y ) =N U * (u x (k) M ;(k)). (3.4) 

According to Eq. (12. 7D this quantity is expected to vanish. This is again confirmed nicely 
in Fig.|3pb) to three orders of magnitude. In contrast measurements of R k (u x ,u y ) in an 
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Figure 3: R^UxUy) averaged over 8xl0 6 iterations in a T fl = l for all a, V = 2Q 2 fluctuating D2Q9 ideal gas 
simulation with and without active ghost noise. Again, take note of the different scales. 

implementation without ghost noise exhibits significant correlations of up to 0.25/5 for 
intermediate k x and k y ranges as seen in Fig.[3la). 

The required condition in Eq. §1.7) , u<l, suggests that this noise implementation 
may suffer from a lack of Galilean invariance. To estimate the magnitude of this violation 
we consider an imposed mean velocity in the x-direction. We measured correlators for 
a fluctuating system with large superimposed velocity of u x = 0.1. The results in Fig. [4] 
indicate that indeed the moment fluctuations do not decouple and Eq. (12.711 is no longer 
fulfilled. Compared to the Ladd implementation these errors are still smaller, but they 
approach the same order of magnitude for maximal accessible velocities. A more com- 
prehensive investigation of these effects is subject of a forthcoming publication. 



4 Discussion and Outlook 

We have shown here that the Adhikari approach to use an improved LB method presents 
a promising scheme to simulate fluctuating hydrodynamics. The ability to interprete the 
ghost degrees of freedom as resulting from discrete particle distributions gives us the 
ability to systematically introduce fluctuations. This approach recovers fluctuations not 
only in the hydrodynamic limit but also for much shorter wavelengths. However, this is 
only true in the absence of flow. Since lattice Boltzmann methods are not generally used 
in this regime one may wonder if Galilean invariance violations may not erase much of 
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Figure 4: Correlators Sk(p), S^Ux), and R^(u x ,Uy) averaged over 5 x 10 6 iterations or a T fl =l for all a, V=20 2 
fluctuating D2Q9 ideal gas simulation with a constant velocity of it* = 0.1. 



the improvement achieved by including noise in the ghost modes. This is a subject to 
which we will return in a forthcoming paper. 
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